4  Feux

4.1 Feux au cours de l’année écoulée

Pour cette visualisation, nous avons mobilisé des données provenant du satellite VIIRS (Visible Infrared Imaging Radiometer Suite). Les données téléchargées depuis le service FIRMS de la NASA portent uniquement sur les douze derniers mois. Elles incluent plusieurs variables, parmi lesquelles nous utilisons principalement la variable ‘frp’ (Fire Radiative Power) pour représenter l’intensité des incendies. Cette variable, exprimée en mégawatts (MW), quantifie le taux d’émission radiative des feux actifs, offrant une indication directe de l’intensité du feu. De plus, nous avons défini une zone d’intérêt de 20km autour de l’aire protégée et nous visualisons les données à l’intérieur de cette zone. La carte finale montre ces incendies avec une gradation de couleurs pour refléter l’intensité de chaque incendie.En cliquant sur les points, on peut consulter la date à laquelle le foyer a été enregistré par le satellite, ainsi que la probabilité qu’il s’agisse téellement d’un feu (classé en “élevée”, “nominale” et “faible”).

Code
# Ci-dessous les librairies R utilisées pour cette analyse
library(tidyverse) # pour la manipulation et visualisation de données
library(sf) # pour traiter les données spatiales
library(tmap) # pour générer des cartes
library(aws.s3) # Pour accéder aux données volumineuses stockées en ligne
library(geodata) # Pour avoir les frontières administratives
library(curl) # Pour télécharge
library(lubridate) # Pour manipuler des dates
library(httr)
library(mapme.biodiversity)
library(units)
library(arrow)
Code
# Create a dataframe of the coordinates
coordinates <- data.frame(
  location = c("Beroroha", "Beronono", "Tsivoko", "Makaykely"),
  latitude = c("21°40.504’S", "21°21.669’S", "21°17.712’S", "21°28.074’S"),
  longitude = c("45°9.571’E", "45°14.885’E", "45°22.732’E", "45°21.896’E")
)

# Function to parse DMS and convert to decimal degrees
ddm_to_decimal <- function(dms){
  # Extraire les degrés, minutes et secondes
  parts <- regmatches(dms, gregexpr("[0-9.]+", dms))[[1]]
  degree <- as.numeric(parts[1])
  minute <- as.numeric(parts[2]) / 60
  direction <- ifelse(grepl("S|W", dms), -1, 1)
  
  # Converstion au format décimal
  decimal <- direction * (degree + minute)
  
  return(decimal)
}

# Apply function to each coordinate
coordinates$latitude_decimal <- sapply(coordinates$latitude, ddm_to_decimal)
coordinates$longitude_decimal <- sapply(coordinates$longitude, ddm_to_decimal)

# Convertir les coordonnées en objet sf
coordinates_sf <- st_as_sf(coordinates, 
                           coords = c("longitude_decimal", "latitude_decimal"), 
                           crs = 4326)


makay_ap <- st_read("data/Makay_wgs84.geojson", quiet = TRUE) %>%
  rename(Statut = SOURCETHM) %>%
  mutate(Statut = recode(Statut, "zone tampon" = "Zone tampon"))

# Fusion des différentes aires
makay_union <- makay_ap %>%
  st_union() %>%
  st_sf() %>%
  st_make_valid()

# Define the buffer size (in kilometers) for the square expansion
buffer_size_km <- 20

# Create a 20km buffer around Makay PA
makay_buffer <- makay_union %>%
  st_buffer(dist = set_units(buffer_size_km, km)) %>%
  st_as_sf()

# Manual request at https://firms.modaps.eosdis.nasa.gov/download/create.php
# Uses the following coordinates : 44.4,-22.4,46.4,-19.9
# For VIIRS NOAA-20
firms_focus <- read_csv("firms_manual.csv") %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
  filter(lengths(st_within(., makay_buffer)) > 0)


breaks_frp <- c(min(firms_focus$frp), 
                median(firms_focus$frp), 
                max(firms_focus$frp))
breaks_frp_5 <- quantile(firms_focus$frp, probs = seq(0, 1, by = 0.25), 
                         na.rm = TRUE)
tmap_mode("view")
tm_shape(firms_focus) +
  tm_dots(col = "frp", size = 1, alpha = 0.2,
          palette = c("blue", "cyan", "yellow", "orange", "red"),
          border.col = NULL, breaks = breaks_frp_5,
          popup.vars = c("acq_date", "confidence", "frp"),
          title = "Fire radiative power") +
  tm_layout(title = "Feux détectés autour du Malay") +
  tm_shape(makay_union) +
  tm_borders(col = "darkblue") + 
  tm_shape(coordinates_sf) +
  tm_dots(col = "black", size = 0.1, labels = coordinates_sf$location) +
  tm_view(set.view = 10) +
  tm_text(text = "location", ymod = -0.5) +
  tm_basemap("OpenStreetMap") +
  tm_scale_bar()

4.2 Feux détectés sur deux décennies

Ces données s’appuient sur les données MODIS qui sont moins précises mais couvrent une période plus longue. On compare le Makay avec une zone périphérique de 60km. This analysis is based on the example produced by Darius Goergen on the Mapme website.

Code
  s3_output <- "diffusion/fire-db2.parquet"
  makay_union <- st_make_valid(makay_union)
if (is.na(st_crs(makay_union))) st_crs(makay_union) <- 4326

# Utilise un CRS métrique pour les buffers (Makay est vers ~45E, 20S -> UTM 38S)
makay_m <- st_transform(makay_union, 32738)    # EPSG:32738 (UTM 38S)
buf20_m  <- st_buffer(makay_m, 20000)          # 20 km
# Donut = buffer - polygone
donut_m  <- st_make_valid(st_difference(buf20_m, st_buffer(makay_m, 0)))

# Reviens en WGS84 pour les jointures avec (x,y)
makay_wgs  <- st_transform(makay_m, 4326)
donut_wgs  <- st_transform(donut_m, 4326)

# 1) Ouvrir le dataset Parquet (Arrow détecte les partitions Hive tout seul)
ds <- open_dataset(paste0("s3://projet-betsaka/", s3_output), format = "parquet")

# 2) Pré-filtrer côté Arrow par bbox + variable (réduction d’E/S)
bb  <- st_bbox(st_union(makay_wgs))
tbl <- ds %>%
  filter(
    variable == "fire_occurrences_1year",
    x >= !!bb["xmin"], x <= !!bb["xmax"],
    y >= !!bb["ymin"], y <= !!bb["ymax"]
    # , year >= 2002  # décommente si tu veux mimer l’exclusion 2000–2001
  ) %>%
  select(pixelid, x, y, year, value) %>%
  collect()

# 3) Conversion en points sf puis tag "zone"
pts <- st_as_sf(tbl, coords = c("x","y"), crs = 4326, remove = FALSE)

in_makay  <- as.logical(st_intersects(pts, makay_wgs, sparse = FALSE)[,1])
in_donut  <- as.logical(st_intersects(pts, donut_wgs, sparse = FALSE)[,1])

pts$zone <- NA_character_
pts$zone[in_makay] <- "Makay"
pts$zone[!in_makay & in_donut] <- "Périphérie 20km"

pts2 <- pts %>% filter(!is.na(zone))

# 4) Agrégats annuels (et total si tu veux)
res_by_year <- pts2 %>%
  st_drop_geometry() %>%
  group_by(zone, year) %>%
  summarise(
    n_cells  = n(),
    sum_fire = sum(value, na.rm = TRUE),
    mean_fire= mean(value, na.rm = TRUE),
    .groups  = "drop"
  ) %>%
  arrange(zone, year)

res_total <- pts2 %>%
  st_drop_geometry() %>%
  group_by(zone) %>%
  summarise(
    n_cells  = n(),
    sum_fire = sum(value, na.rm = TRUE),
    mean_fire= mean(value, na.rm = TRUE),
    .groups  = "drop"
  )

res_by_year
# A tibble: 50 × 5
   zone   year n_cells sum_fire mean_fire
   <chr> <int>   <int>    <dbl>     <dbl>
 1 Makay  2000    3364       70    0.0208
 2 Makay  2001    3364      135    0.0401
 3 Makay  2002    3364      318    0.0945
 4 Makay  2003    3364      416    0.124 
 5 Makay  2004    3364      274    0.0815
 6 Makay  2005    3364      601    0.179 
 7 Makay  2006    3364      337    0.100 
 8 Makay  2007    3364      391    0.116 
 9 Makay  2008    3364      411    0.122 
10 Makay  2009    3364      315    0.0936
# ℹ 40 more rows
Code
res_total
# A tibble: 2 × 4
  zone            n_cells sum_fire mean_fire
  <chr>             <int>    <dbl>     <dbl>
1 Makay             84100     7467    0.0888
2 Périphérie 20km   95175    10754    0.113 

Compute indicators on fires.

Code
p_sum <- ggplot(res_by_year, aes(x = year, y = sum_fire, color = zone)) +
  geom_line() +
  geom_point(size = 1) +
  labs(x = "Année", y = "Somme annuelle de feux",
       color = "Zone",
       title = "Somme annuelle de feux (Makay vs Périphérie 20 km)") +
  theme_minimal()

# 2) Moyenne annuelle par zone (courbe)
p_mean <- ggplot(res_by_year, aes(x = year, y = mean_fire, color = zone)) +
  geom_line() +
  geom_point(size = 1) +
  labs(x = "Année", y = "Moyenne annuelle de feux par cellule",
       color = "Zone",
       title = "Moyenne annuelle de feux par cellule (Makay vs Périphérie 20 km)") +
  theme_minimal()

# 3) Moyenne sur l'ensemble des années (barres)
avg_by_zone <- res_by_year %>%
  group_by(zone) %>%
  summarise(
    mean_sum_fire  = mean(sum_fire, na.rm = TRUE),
    mean_mean_fire = mean(mean_fire, na.rm = TRUE),
    .groups = "drop"
  )

# a) moyenne annuelle de la somme (barres)
p_avg_sum <- ggplot(avg_by_zone, aes(x = zone, y = mean_sum_fire, fill = zone)) +
  geom_col() +
  labs(x = NULL, y = "Somme annuelle moyenne",
       title = "Somme annuelle moyenne de feux (2000–2024)") +
  theme_minimal() +
  theme(legend.position = "none")

# b) moyenne (sur années) de la moyenne par cellule (barres) — optionnel
p_avg_mean <- ggplot(avg_by_zone, aes(x = zone, y = mean_mean_fire, fill = zone)) +
  geom_col() +
  labs(x = NULL, y = "Moyenne annuelle par cellule",
       title = "Moyenne annuelle par cellule (2000–2024)") +
  theme_minimal() +
  theme(legend.position = "none")

# Afficher
p_sum

Fire occurence count in Makay PA (source: MODIS)
Code
p_mean

Fire occurence count in Makay PA (source: MODIS)
Code
p_avg_sum

Fire occurence count in Makay PA (source: MODIS)

Les différents ordres de grandeur doivent être interprétés avec précaution car la zone périphérique est plus de deux fois plus vaste que l’aire protégée (3,751 km2 contre 8,249 km2).